function [PE_out] = Fn_EstimatePE_Mar2022(A_hat, W_ij, y1, Cons)
    %alpha_educ, alpha_gender, Ehat_E1_educ, Ehat_E1_gender, beta_educ, beta_gender Ehat_E2_educ, Ehat_E2_gender] = 
%% Input: W_ij, A_hat
%% Output: Estimates and residuals for PE parameters


% Setup
LMH_educ = [Cons.L_educ, Cons.M_educ, Cons.H_educ];
X_0 = [ones(Cons.n,1), Cons.P, W_ij*Cons.P];
X_0_educ_LMH = [];
for z=1:3
    X_0_educ_LMH = [X_0_educ_LMH, bsxfun(@times,LMH_educ, X_0(:,z))];
end

LMH_gender = [Cons.L_gender, Cons.M_gender, Cons.H_gender];
X_0 = [ones(Cons.n,1), Cons.P, W_ij*Cons.P];
X_0_gender_LMH = [];
for z=1:3
    X_0_gender_LMH = [X_0_gender_LMH, bsxfun(@times,LMH_gender, X_0(:,z))];
end

%%%%% Without A's
    X_educ_m1 = X_0_educ_LMH;
    PE_out.alpha_educ_m1 = inv(X_educ_m1'*X_educ_m1)*X_educ_m1'*y1.y1_educ_m1;
    PE_out.Ehat_educ_m1 = y1.y1_educ_m1 - X_educ_m1*PE_out.alpha_educ_m1;
    PE_out.Sigma2_educ_m1 = PE_out.Ehat_educ_m1'*PE_out.Ehat_educ_m1/Cons.n;

    X_educ_m3 = [X_0_educ_LMH, bsxfun(@times,LMH_educ,Cons.y0_educ), bsxfun(@times,LMH_educ,W_ij*Cons.y0_educ)];
    PE_out.alpha_educ_m3 = inv(X_educ_m3'*X_educ_m3)*X_educ_m3'*y1.y1_educ_m3;
    PE_out.Ehat_educ_m3 = y1.y1_educ_m3 - X_educ_m3*PE_out.alpha_educ_m3;
    PE_out.Sigma2_educ_m3 = PE_out.Ehat_educ_m3'*PE_out.Ehat_educ_m3/Cons.n;


    X_gender_m1 = X_0_gender_LMH;
    PE_out.alpha_gender_m1 = inv(X_gender_m1'*X_gender_m1)*X_gender_m1'*y1.y1_gender_m1;
    PE_out.Ehat_gender_m1 = y1.y1_gender_m1 - X_gender_m1*PE_out.alpha_gender_m1;
    PE_out.Sigma2_gender_m1 = PE_out.Ehat_gender_m1'*PE_out.Ehat_gender_m1/Cons.n;

    X_gender_m3 = [X_0_gender_LMH, bsxfun(@times,LMH_gender,Cons.y0_gender), bsxfun(@times,LMH_gender,W_ij*Cons.y0_gender)];
    PE_out.alpha_gender_m3 = inv(X_gender_m3'*X_gender_m3)*X_gender_m3'*y1.y1_gender_m3;
    PE_out.Ehat_gender_m3 = y1.y1_gender_m3 - X_gender_m3*PE_out.alpha_gender_m3;
    PE_out.Sigma2_gender_m3 = PE_out.Ehat_gender_m3'*PE_out.Ehat_gender_m3/Cons.n;
 



%%%%% WITH A'S
if var(A_hat)>1e-4
    X_educ_m2 = [X_0_educ_LMH, bsxfun(@times,LMH_educ,A_hat), bsxfun(@times,LMH_educ,W_ij*A_hat)];
    PE_out.alpha_educ_m2 = inv(X_educ_m2'*X_educ_m2)*X_educ_m2'*y1.y1_educ_m2;
    PE_out.Ehat_educ_m2 = y1.y1_educ_m2 - X_educ_m2*PE_out.alpha_educ_m2;
    PE_out.Sigma2_educ_m2 = PE_out.Ehat_educ_m2'*PE_out.Ehat_educ_m2/Cons.n;

    X_educ_m4 = [X_0_educ_LMH, bsxfun(@times,LMH_educ,A_hat), bsxfun(@times,LMH_educ,W_ij*A_hat), bsxfun(@times,LMH_educ,Cons.y0_educ), bsxfun(@times,LMH_educ,W_ij*Cons.y0_educ)];
    PE_out.alpha_educ_m4 = inv(X_educ_m4'*X_educ_m4)*X_educ_m4'*y1.y1_educ_m4;
    PE_out.Ehat_educ_m4 = y1.y1_educ_m4 - X_educ_m4*PE_out.alpha_educ_m4;
    PE_out.Sigma2_educ_m4 = PE_out.Ehat_educ_m4'*PE_out.Ehat_educ_m4/Cons.n;
    
    
    X_gender_m2 = [X_0_gender_LMH, bsxfun(@times,LMH_gender,A_hat), bsxfun(@times,LMH_gender,W_ij*A_hat)];
    PE_out.alpha_gender_m2 = inv(X_gender_m2'*X_gender_m2)*X_gender_m2'*y1.y1_gender_m2;
    PE_out.Ehat_gender_m2 = y1.y1_gender_m2 - X_gender_m2*PE_out.alpha_gender_m2;
    PE_out.Sigma2_gender_m2 = PE_out.Ehat_gender_m2'*PE_out.Ehat_gender_m2/Cons.n;

    X_gender_m4 = [X_0_gender_LMH, bsxfun(@times,LMH_gender,A_hat), bsxfun(@times,LMH_gender,W_ij*A_hat), bsxfun(@times,LMH_gender,Cons.y0_gender), bsxfun(@times,LMH_gender,W_ij*Cons.y0_gender)];
    PE_out.alpha_gender_m4 = inv(X_gender_m4'*X_gender_m4)*X_gender_m4'*y1.y1_gender_m4;
    PE_out.Ehat_gender_m4 = y1.y1_gender_m4 - X_gender_m4*PE_out.alpha_gender_m4;
    PE_out.Sigma2_gender_m4 = PE_out.Ehat_gender_m4'*PE_out.Ehat_gender_m4/Cons.n;    


end